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Abstract 

Some of the most interesting Higgs-production processes at future e + e~ 
colliders are of the type e + e~ — > ffH. We present a calculation of the 
complete O(a) corrections to these processes in the Standard Model for 
final-state neutrinos and top quarks. Initial-state radiation beyond O(a) at 
the leading-logarithmic level as well as QCD corrections are also included. 
The electroweak corrections turn out to be sizable and reach the order of 
±10% and will thus be an important part of precise theoretical predictions 
for future e + e~ colliders. 

Furthermore, an overview is given of a technique for a fast and reliable 
numerical calculation of multi-leg one- loop integrals. The method is numer- 
ically stable also for exceptional momentum configurations and easily allows 
the introduction of complex masses and the calculation of higher orders in 
the expansion around D = 4. 

1 Introduction 

One of the main future tasks in particle physics will be the investigation of the 
mechanism of electroweak symmetry breaking in general and the discovery of the 
Higgs boson and the determination of its properties in particular. Since the Higgs- 
boson mass is expected to be in the range from the lower experimental bound of 
114.4 GeV up to 1 TeV, with a light Higgs mass (below ~ 200 GeV) favoured by 
electroweak precision data, the LHC will be able to discover it in the full mass 
range, provided it exists and has no exotic properties. However, for the complete 
determination of its profile, including its couplings to fermions and gauge bosons, 
experiments in the clean environment of an e + e~ linear collider are indispensable. 

Here we concentrate on the associated production of a Higgs boson together 
with a pair of neutrinos or top quarks in e + e~ annihilation, which are among the 
most interesting Higgs-boson production processes at future e + e~ linear colliders. 

*Talk given at the final meeting of the European Network "Physics at Colliders", Montpellier, 
September 26-27, 2004. 
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The calculation of the radiative corrections to these processes is presented in the 
next two sections. 

The last section gives a sketch of a technique for a fast and reliable numerical 
calculation of multi-leg one-loop integrals and describes an implementation of the 
method in Mathematica and C++. 

2 The process e + e _ — > vvtl 

At e + e~ colliders the two main Higgs production processes are the Higgs-strahlung 
and W-boson-fusion processes. In the Higgs-strahlung process the Higgs boson is 
radiated off a Z boson, with the corresponding cross section rising sharply at 
the threshold, located at a centre-of-mass (CM) energy of y/s = M z + M H , to a 
maximum a few tens of GeV above the threshold energy and then falling off as 
1/s. In the W-boson-fusion process the Higgs boson is produced via fusion of two 
W bosons, each emitted from an incoming electron/positron. The corresponding 
cross section grows as In s and thus is the dominant production mechanism at large 
energies. Both production mechanisms appear in the process e + e _ — > z^H, with 
/ = e, fi, or r, though the W-boson-fusion process is only present for I = e. 

For the process e + e~ — > ZH the 0(a) electroweak radiative corrections have 
been calculated many years ago in Ref. [1]. Furthermore a Monte Carlo algorithm 
for the calculation of the real photonic corrections to this process was described 
in Ref. [2]. For the full process e + e~ — ► z/z/H there has been a lot of activity re- 
garding the electroweak corrections recently. Within the Minimal Supersymmetric 
Standard Model (MSSM) the fermion and sfermion loop contributions have been 
evaluated in Refs. [3, 4]. Analytical results for the one-loop corrections in the SM 
have been obtained in Ref. [5], though no numerical results have been given there. 
Finally, calculations of the complete 0(a) electroweak corrections to e + e" — > z/z/H 
in the SM have been performed in Refs. [6, 7]. Very recently also results on correc- 
tions to the Z-boson-fusion process e + e~ — > e + e~H have been presented in Ref. [8]. 

2.1 Calculational framework 

The calculation of the one-loop diagrams has been carried out in the 't Hooft- 
Feynman gauge using standard techniques. The renormalization is carried out in 
the on-shell renormalization scheme, as e.g. described in Ref. [9]. The electron 
mass m c is neglected whenever possible. 

The calculation of the Feynman diagrams has been performed in two com- 
pletely independent ways, leading to two independent computer codes for the 
numerical evaluation. Both calculations are based on the methods described in 
Ref. [9]. Apart from the 5-point functions the tensor coefficients of the one-loop 
integrals are recursively reduced to scalar integrals with the Passarino-Veltman 
algorithm [10] at the numerical level. The scalar integrals are evaluated using the 
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methods and results of Refs. [9, 11], where ultraviolet divergences are regulated 
dimensionally and IR divergences with an infinitesimal photon mass. The 5-point 
functions are reduced to 4-point functions following Ref. [12], where a method for a 
direct reduction is described that avoids leading inverse Gram determinants which 
potentially cause numerical instabilities. As a check of gauge independence the 
calculation of the virtual corrections has been repeated using the background-field 
method [13]. 

The results of the two different codes, and also those obtained within the 
conventional and background-field formalism, are in good numerical agreement 
(typically within at least 12 digits for non-exceptional phase-space points). 

We use two different schemes for the inclusion of the finite Z-boson decay 
width. In the fixed-width scheme, each resonant Z-boson propagator l/(s uP — M|), 
where s vV is the invariant mass of the neutrino-antineutrino pair, is replaced by 
l/Cw — M% + iM z T z ), while non-resonant contributions are kept untouched. This 
potentially violates gauge invariance, because the resonant part of the amplitude 
alone is not gauge invariant. As a second option, we applied a factorization scheme 
where the full (gauge-invariant) ZH-production amplitude with zero Z-boson width 
is rescaled by a factor (s uP — M|) / (s uP — M| + iMzTz). However in this scheme the 
non-resonant part of the ZH-amplitude is neglected on resonance. Nevertheless 
both schemes give the same results for the total cross section within integration 
errors. 

The matrix elements for the real photonic corrections are evaluated using the 
Weyl-van der Waerden spinor technique as formulated in Ref. [14] and have been 
successfully checked against the result obtained with the package Madgraph [15]. 
The soft and collinear singularities are treated both in the dipole subtraction 
method following Refs. [16, 17] and in the phase-space slicing method following 
closely Ref. [18]. 

The emission of photons collinear to the incoming electrons or positrons leads 
to corrections that are enhanced by large logarithms of the form \n(m 2 c /s). In 
order to achieve an accuracy at the few 0.1% level, the corresponding higher-order 
contributions, i.e. contributions beyond O(o), must be taken into account. These 
are included in our calculation at the leading-logarithmic level using the structure 
functions given in Ref. [19] (for the original papers see references therein). 

The calculation is done in the so-called G^-scheme, i.e. we derive the elec- 
tromagnetic coupling a = e 2 / (Air) from the Fermi constant according to 
q;g m = V^G^M^s^/tt. This procedure absorbs the corrections proportional to 
m 2 /M^j in the fermion-W-boson couplings and the running of ct(Q 2 ) from Q 2 = 
to the electroweak scale. In the relative radiative corrections, we use a(0) as cou- 
pling parameter, which is the correct effective coupling for real photon emission. 

The cross section for e + e~ — > uuYL is dominated by the WW-fusion diagram, 
which gets its main contribution from the region of small momentum transfers. 
Consequently, the corresponding corrections are determined by the ez/ e W and 
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WWH vertex corrections for small invariant W masses. The correction to the 
ez/ e W vertex and the main contributions to the WWH vertex in the relevant kine- 
matical region are well approximated by Ar. Thus, parametrizing the lowest 
order in terms of (G M -scheme) absorbs a large part of the universal correc- 
tions. Further universal corrections have been obtained by extracting the leading 
mt-dependent corrections of the WW-contribution in the heavy-top limit in the 
G M -scheme. These reproduce the full m t -dependent corrections rather well for the 
WW channel, which is dominated by small momentum transfers. Therefore, we 
have defined the following improved Born approximation (IBA) 

i non-photonic 1 i WW m t /i\ 

iba =da -da ^JfT- (1) 

The corresponding expression for the m t — > oo limit of the ZH contribution is not 
included in the definition of the IBA, since it does not give a good description. In 
the ZH channel y^i is a typical scale for the momentum transfer, which is larger 
than m t in the physically interesting region of e + e~ — > z/z/H. Finally, daj 1 ^" 15110 * 01110 
is convoluted with the ISR structure functions to yield the cross section of the full 
IBA. 

The phase-space integration is performed with Monte Carlo techniques in both 
computer codes. The first code employs a multi-channel Monte Carlo generator 
similar to the one implemented in RacoonWW [17, 20] and Lusifer [21], the second 
one uses the adaptive multi-dimensional integration program Vegas [22]. 



2.2 Comparison to related work 

We have compared our results for the O{o) corrections to Ref. [7] and the contri- 
butions from closed fermion loops with Refs. [3, 4]. 

Adapting the input parameters and the parametrization of the lowest-order 
matrix element to those used by Belanger et al. [7], we reproduced the numbers 
for the total cross section given in Table 2 of the first paper of Ref. [7]. Note that 
we switch off the ISR beyond 0(a) in this comparison. In Table 1 we list for each 
Higgs-boson mass the results of Ref. [7] 1 together with our results. The numbers 
in parenthesis indicate the errors in the last digits. We find agreement within 10~ 4 
for the total lowest-order cross section and within 0.3% for the corrected cross 
section. The corrections relative to the lowest-order cross section agree within 
0.2%. This is of the order of the statistical error of Ref. [7], which is about 0.1%. 
Note that Belanger et al. use a(0) to parametrize the lowest-order cross section. 
As a consequence their relative corrections are shifted by 3Ar ss +9% compared 
to those in the G^-scheme. 

1 According to F. Boudjcma, the numbers for the lowest-order cross section in Table 2 of 
Ref. [7] have integration errors of the order of 0.2%. Table 1 contains updated numbers obtained 
with increased statistics. 
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M H [GeV] 


0"tree [ft>] 


<7[fb] 


6[%] 




150 


61.074(7) 


60.99(7) 


-0.2 


Ref. [7] 




61.076(5) 


60.80(2) 


-0.44(3) 


this work 


250 


21.135(2) 


20.63(2) 


-2.5 


Ref. [7] 




21.134(1) 


20.60(1) 


-2.53(3) 


this work 


350 


4.6079(5) 


4.184(4) 


-9.1 


Ref. [7] 




4.6077(2) 


4.181(1) 


-9.27(3) 


this work 



Table 1: Total cross section in lowest order and including the full 0(a) corrections 
and the relative corrections for yfs = 500 GeV and various Higgs masses for the 
input parameter scheme of Ref. [7] 

We have also reproduced the cos#h and En distributions in Figures 1 and 2 
of the first paper of Ref. [7]. We found agreement within the accuracy of these 
figures. 

When considering only fermion-loop corrections, we find agreement with the 
calculations of Refs. [3, 4] , once the appropriate renormalization and input-parame- 
ter schemes are adopted. For more details on this comparison we refer to Ref. [6]. 

2.3 Numerical results 

The results for the total cross section in lowest order and including the radiative 
corrections are shown in Figure 1 on the l.h.s. as a function of the CM energy for 
Mh = 150 GeV. The relative corrections shown on the r.h.s are large (< —20%) 
and vary strongly in the ZH-threshold region while they are flat and about —10% 
for energies above 500 GeV. They are always negative because they are dominated 
by initial-state radiation and the cross section is monotonously rising. Also shown 
in Figure 1 on the r.h.s are the residual relative corrections normalized to the IBA 
which are about 1% near the threshold and reach 3-4% at high energies. Although 
they are systematically smaller than the corrections relative to the lowest order 
in the scheme, the inclusion of the full 0(a) corrections is necessary for a 
precision analysis. 

3 The process e + e _ — > ttH 

We have also investigated the process e + e~ — > ttH, which is interesting since it 
permits a direct access to the top-quark Yukawa coupling #ttH> which is by far 
the largest Yukawa coupling (g t m ~ 0.5) in the SM. This is possible because 
the process proceeds mainly through Higgs-boson emission off top quarks, while 
emission from intermediate Z bosons plays only a minor role if the Higgs-boson 
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Figure 1: Lowest-order and corrected cross sections (l.h.s.) as well as relative 
corrections with respect to Born result and improved Born approximation (r.h.s.) 
in the G M scheme for a Higgs-boson mass Mh = 150 GeV 



mass is not too large, i.e. Mh ~ 100-200 GeV. For a light Higgs boson with a 
mass around Mh ~ 120 GeV, a precision of about 5% can be reached at an e + e~ 
linear collider operating at y/s = 800 GeV with a luminosity of J L dt ~ 1000 fb _1 
[23] . An even better accuracy can be obtained by combining the ttH channel with 
information from other Higgs-production and decay processes in a combined fit 
[24]- 

Within the SM the 0(a s ) corrections have been calculated for the dominant 
photon-exchange channel in Ref. [25], while the full set of diagrams has been 
evaluated in Ref. [26]. The 0(a s ) corrections to the photon-exchange channel in 
the MSSM have been considered in Ref. [27]. In Ref. [28] all QCD diagrams have 
been taken into account, while the SUSY-QCD corrections have been worked out in 
Ref. [29]. The evaluation of the electroweak 0(a) corrections in the SM has made 
considerable progress recently. Results have been presented in Refs. [30, 31, 32], 
with agreement between Refs. [31, 32] while Ref. [30] shows deviations close to 
threshold and at high energies. 

Our calculation [32] includes the 0(a) electroweak and the 0(a s ) QCD cor- 
rections. Though the calculation of the virtual corrections for this process is much 
more involved than for the process e + e~ — > i/z/H, the same calculational techniques 
could be used. 
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v^[GeV] 


Otree [ft>] 


<j[fb] 


*[%] 




500 


4.8142- 10~ 4 


3.401 • 10" 4 


-29.35 


Ref. [30] 




4.8140(8) • 10~ 4 


3.168(4) • 10- 4 


-34.19(8) 


this work 


800 


1.58 


1.63 


3.60 


Ref. [30] 




1.5749(2) 


1.6243(4) 


3.14(2) 


this work 


1000 


1.47 


1.53 


4.47 


Ref. [30] 




1.4664(2) 


1.5273(4) 


4.15(2) 


this work 


2000 


0.6270 


0.6297 


0.43 


Ref. [30] 




0.6269(1) 


0.6526(3) 


4.11(5) 


this work 



Table 2: Total cross section in lowest order and including the full electroweak 
O(a) corrections as well as the relative corrections for Mh = 150 GeV and various 
CM energies for the input-parameter scheme of Ref. [30] . The statistical errors of 
Ref. [30] are estimated by the authors to be below 1% 

3.1 Comparison to related work 

The results on the QCD corrections have been reproduced with the (publically 
available) computer code based on the calculation of Ref. [26] . We found agreement 
within the statistical integration errors. 

For a comparison of the electroweak 0(a) corrections with the results of Ref. [30] 
we changed our input parameters to the ones quoted there and switched to the 
a(0)-scheme. In Table 2 we compare some representative numbers 2 from the calcu- 
lation of Ref. [30] with the corresponding results from our Monte Carlo generator. 
The numbers in parentheses give the errors in the last digits of our calculation. 
The tree-level cross sections coincide within 0.03%. Most of the numbers for the 
one-loop corrected cross sections agree within 1-2%, i.e. roughly within the esti- 
mated error of Ref. [30] . However, for the corrected cross sections at y/s = 2 TeV, 
i.e. at high energies, and the one very close to threshold, i.e. for y/s = 500 GeV 
and Mh = 150 GeV, we find differences of 4% and 7%, respectively. The same 
holds for the relative corrections. Ours are larger by about 4% at \/s = 2 TeV and 
smaller by about 5% for the selected cross section close to threshold. 

Finally, we have also compared the electroweak 0{a) corrections with Ref. [31], 
where the a(0)-scheme has been used. In Table 3 we list the results of Table 2 
of Ref. [31] for Mh = 120 GeV together with the corresponding results from our 
Monte Carlo generator. Again the numbers in parentheses give the errors in the 
last digits. We reproduce the results for the lowest-order cross section within the 

2 These numbers were kindly provided to us by Zhang Ren- You and You Yu quoting a statis- 
tical error below 1%. 
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Vi[GeV] 


Otree [ft>] 


<j[fb] 


5[%] 




600 


1.7293(3) 


1.738(2) 


0.5 


Ref. [31] 




1.7292(2) 


1.7368(6) 


0.44(3) 


this work 


800 


2.2724(5) 


2.362(4) 


3.9 


Ref. [31] 




2.2723(3) 


2.3599(6) 


3.86(2) 


this work 


1000 


1.9273(5) 


2.027(4) 


5.2 


Ref. [31] 




1.9271(3) 


2.0252(5) 


5.09(2) 


this work 



Table 3: Total cross section in lowest order and including the full electroweak 0(a) 
corrections as well as the relative corrections for Mh = 120 GeV and various CM 
energies for the input-parameter scheme of Ref. [31]. 

integration errors, which are about 2-3 x 10~ 4 . The results for the cross section 
including electroweak corrections as well as the relative corrections coincide to 
better than 0.1% which is of the order of the integration error of the results of 
Ref. [31]. 

3.2 Numerical results 

Results for the total cross section in lowest order and the corrected cross section 
including both the electroweak and QCD corrections are shown in Figure 2 on the 
l.h.s. Away from the kinematic threshold at y/s = 2m t + M H the size of the cross 
section is typically a few fb, with a maximum at about 800 GeV. On the r.h.s. 
of Figure 2 the relative corrections are shown. The QCD corrections are large 
and positive close to threshold where soft-gluon exchange in the tt system leads 
to a Coulomb-like singularity. For larger energies the QCD corrections decrease, 
eventually turn negative and reach about —8% at an energy of y/s = 1.5 TeV. 
The electroweak corrections are about —10% and vary only weakly with energy 
away from the threshold region, and are thus of a comparable size as the QCD 
corrections. Close to threshold they reach about —20% due to the large ISR QED 
corrections in this region. The behaviour of the combined electroweak and QCD 
corrections is dominated by the Coulomb-like singularity close to threshold while 
turning negative and reaching about —15% at high energies. 

Summarizing, for both of the processes e + e~ — > vvH and e + e~ — > ttH the 
0(a) corrections are sizeable and typically of the order ±10%. They will thus 
be an important ingredient of precise theoretical predictions for future e + e~ col- 
liders. Our results agree with the ones of an independent calculation within the 
integration errors, which are around 0.1-0.2%. Moreover, these calculations show 
that techniques for the calculation of one-loop corrections to 2 — > 3 processes are 
available and work well in practical applications. 
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Figure 2: Lowest-order and corrected cross sections (l.h.s.) as well as relative 
corrections (r.h.s.) in the scheme for a Higgs-boson mass Mh = 150 GeV 



4 Numerical calculation of one-loop integrals 

In this section we present a technique for a fast and reliable numerical calcu- 
lation of multi-leg one-loop integrals and describe an implementation in Mathe- 
matica/C-|— h The method is numerically stable also for exceptional momentum 
configurations and easily allows the introduction of complex masses and the cal- 
culation of higher orders in the expansion around D = 4. 

Using the conventional analytic approach of Ref. [11] all scalar loop integrals 
can be expressed in terms of dilogarithms and logarithms. Furthermore, using the 
reduction algorithm of Ref. [10] all tensor loop integrals, i.e. integrals containing 
loop momenta in the numerator, can be expressed in terms of scalar integrals. 
Therefore, a full analytic solution for one-loop integrals exists. However this ap- 
proach has a number of drawbacks. First of all, with an increasing number of 
external legs the number of dilogarithms in the analytic expression of a scalar 
integral increases rapidly. This can lead to cancellations for multi-leg integrals in 
certain kinematic regions [33]. Furthermore the tensor reduction of Ref. [10] intro- 
duces inverse Gram determinants. These can vanish at the phase space boundary 
even though the tensor coefficients themselves are regular in this region. There 
are thus cancellations among terms in the numerator that can lead to numerical 
instabilities. Unstable particles are also an important issue in multi-leg loop calcu- 
lations, since they appear as virtual particles in the diagrams. One way of dealing 
with them is the introduction of complex masses for the unstable particles. This 
requires the evaluation of loop integrals with complex masses which is cumbersome 
in analytic calculations. Finally, within dimensional regularization the evaluation 



9 



of the loop-by-loop contribution to a 2-loop calculation makes it necessary to ex- 
pand the one-loop integrals beyond the constant term in the expansion around 
D = 4. An analytic calculation of these higher-order terms is rather complicated 
It seems therefore worthwhile to explore alternative numerical approaches to 
the evaluation of one-loop tensor and scalar integrals. The strategy adopted here 
is described in detail in Ref. [34]. It is based on the Bernstein-Tkachov theorem 
[35] which can be used to rewrite one-loop integrals in Feynman-parametric rep- 
resentation in a form better suited for numerical evaluation. The general method 
is outlined in the next section and a description of an implementation in Mathe- 
matica and C++ is given in the last section. 

4.1 Description of the method 

Within dimensional regularization in D = 4 — 2e dimensions any scalar one-loop 
integral can be expressed as an integration over Feynman parameters 

jo = (^r D f d n q 1 

N I?!" 2 J [q 2 - mj][(q + pi) 2 - m 2 .] ■ ■ ■ [(q + p N ~ 1 ) 2 - m 2 N ] 

= (4vr/i 2 r T(N - 2 + e)(-l) N J dS N ^V(x t r (N - 2+e) (2) 
where the integration over Feynman parameters is defined as 

dS n = dxi dx 2 ■ ■■ dx n 
Jo Jo Jo 

and V is a quadratic form in the N — 1 Feynman parameters Xj 

V(x) = x T Hx + 2K T x + L — i5. 

The coefficients H, K and L of V are given in terms of the momenta Pi and the 
masses m«. Note that we use dimensional regularization not only for ultraviolet 
but also for infrared (IR) and collinear singularities. 

In general the quadratic form V can vanish within the integration region, 
though the zero's are shifted into the complex plane by the small imaginary part 
i5. Since the limit 5 — > has to be taken in the end, the form given above is not 
suited for a direct numerical integration. 

Instead, the integral can be rewritten before attempting a numerical evalua- 
tion using the Bernstein-Tkachov theorem [35]. Applied to the case of one-loop 
integrals it states that for any quadratic form V(x) raised to any real power f3 



_ (x - X)jdj 
2(1 + (3) 



V 1+ %x t ) = B-V^ Xl ), (3) 
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where X = —K T H 1 , B = L — K T H 1 K and di = d/dxi. Inserting this relation 
into a Feynman-parameter integral and integrating by parts one obtains 



/ dS n Vf> = 2B{ l +(3) (2 + n + 2(3) J dS n V^ - J dS n ^ g Xl V^ 
where Xi — ~ with X — 1 and X n+1 = and 



(4) 



V(l,xi, . . . ,x n -i) for i = 

Vi(x 1: x n -i) = <( V(xi, . . . , Xi, Xi, . . . , x n -i) for < % < n 
V(xi, . . . , x n _i, 0) for % = n 

Applied to the one-loop integral (2) the first term inside brackets in (4) corresponds 
to the iV-point integral in D + 2 dimensions, while the last term is a sum over 
(N — l)-point integrals in D dimensions obtained by pinching one propagator. 

Recursive application of (4) allows to express any scalar one-loop integral 
as a linear combination of terms of the form J dSk V(xi) m ~ € with any integer 
m > 0. A Taylor expansion up to 0(e a ) will then result in terms of the form 
J dSk V m ■ log 1+a V. For m = the integrand still contains an integrable (loga- 
rithmic) singularity while it is smooth for m > 0. Although larger values of m will 
lead to smoother integrands, the expressions also grow larger due to the repeated 
application of the BT identity (4). The optimal choice for m depends on the chosen 
numerical integration routine and its ability to deal with integrable singularities. 
Note that the calculation of higher orders of the e expansion is straightforward in 
this approach. Furthermore complex masses can also be introduced easily. 

If the integral is infrared or collinear divergent, the repeated application of 
the BT-identity (4) will eventually result in divergent 3-point integrals. For these 
B = and using a modified identity the singularities are automatically extracted 
as poles in e. 

In the case of tensor integrals the parametric representation of the integral 
contains in general Feynman parameters in the numerator. The procedure outlined 
above can also be applied in this case so that no separate reduction to scalar 
integrals is needed. Furthermore, no inverse Gram determinants are introduced 
using this approach, making it numerically reliable also for exceptional kinematic 
configurations. 



4.2 Implementation 

The method outlined above has been implemented in Mathematica and C++ 
with an emphasis on the full automatization of the whole procedure. The user 
only has to supply the algebraic values of the Lorentz invariants calculated from 
the external momenta and the internal masses of the integral. As a result a set of 
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Figure 3: Box-diagram for heavy quark pair production. 



C++ routines with a simple interface is generated. These can then be used for a 
numerical evaluation. 

The implementation first generates the parametric representation for the ten- 
sor coefficients for a given integral up to the maximum desired tensor rank. The 
tensor coefficients are defined according to the conventions of Ref. [9]. In the next 
step consecutive applications of the BT-identities raise the powers of the quadratic 
forms. IR and collinear singularities show up as poles in e during this procedure. 
Then the expansion in e is performed. The power in e up to which the integrals 
are expanded can be chosen by the user. The results of this last algebraic step 
are a number of Feynman-parameter integrals of different dimensions and in some 
cases additional constant terms. Each of the integrands is a vector with compo- 
nents corresponding to the tensor coefficients and the components themselves are 
truncated power series in e. 

The last step is the generation of C++ routines for the calculation of the 
various integrands. Furthermore, a driver routine is generated that performs the 
necessary initializations, calls the numerical integration code and constructs the 
results for the tensor coefficients from the results of the numerical integrations. 
This driver routine is the only part of the code the user interacts with directly. It 
needs only the values for the Lorentz invariants and masses as input and returns 
the coefficients of the e expansion of the tensor coefficients. 

As an example we consider the integral shown in Figure 3, which is both IR and 
collinear divergent. It has been evaluated up to the constant term in the context of 
the calculation of the next-to-leading order QCD corrections to heavy quark pair 
production at hadron colliders [36]. Recently also the 0(e) coefficient has been 
calculated analytically in Ref. [37]. Using our numerical program we obtain for 
s = -t = (500 GeV) 2 and m = 175 GeV 

A)= e- 2 -(-2.85078- 10- 11 +i • 0) 



+ e- 1 ^ 3.87554(7) • lO" 10 - 
+ 1 -(-2.49772(5) • 10~ 9 + 
+ e 1 •( 9.9934(2) • 10~ 9 - 



•4.47800- 10~ n ) 
•6.6096(3) • 10- 10 ) 
•4.4041(2) • 10~ 9 ) 
+ e 2 -(-2.78060(4) • 10~ 8 +i • 1.81859(4) • 10~ 8 ) 
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£W = 9v [ 1 -(-4.59(4) • 10- 7 -i • 5.019(5) • 1(T 6 ) 
+ e 1 -(+1.228(2) • l(T 5 +i • 4.810(2) • lO" 8 ) 
+ e 2 -(-6.697(6) • lO" 5 — i • 2.259(1) • 10" 4 )] 

+ . . . 

where numerical integration errors in the last digit are given in parentheses. These 
results agree with the analytical results within integration errors. 

Our implementation is currently capable of handling all triangle integrals up to 
tensor rank 3 and all box integrals up to rank 4 including IR and collinear divergent 
integrals. All of these can be calculated up to 0(e 2 ). A comparison of the finite part 
and the IR pole with the results of the Loop Tools integral library [38] has shown 
numerical agreement within integration errors for all tensor coefficients of the 3- 
and 4-point functions. For the 5- and 6-point functions only the scalar integrals 
are available so far. The implementation of the remaining tensor coefficients is 
expected to be finished in the near future. 
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